Title: Cytometry-based single cell analysis of intact epithelial signaling reveals MAPK activation divergent from TNF-α-induced apoptosis in vivo

Understanding heterogeneous cellular behaviors in a complex tissue requires the evaluation of signaling networks at single cell resolution. However, probing signaling in epithelial tissues using cytometry-based single cell analysis has been confounded by the necessity of single cell dissociation, where disrupting cell-to-cell connections inherently perturbs native cell signaling states. Here, we demonstrate a novel strategy (Disaggregation for Intracellular Signaling in Single Epithelial Cells from Tissue – DISSECT) that preserves native signaling for Cytometry Time-of-Flight (CyTOF) and fluorescent flow cytometry applications. A 21-plex CyTOF analysis encompassing core signaling and cell-identity markers was performed on the small intestinal epithelium after systemic tumor necrosis factor-alpha (TNF-α) stimulation. Unsupervised and supervised analyses robustly selected signaling features that identify a unique subset of epithelial cells that are sensitized to TNF-α-induced apoptosis in the seemingly homogeneous enterocyte population. Specifically, p-ERK and apoptosis are divergently regulated in neighboring enterocytes within the epithelium, suggesting a mechanism of contact-dependent survival. Our novel single cell approach can broadly be applied, using both CyTOF and multiparameter flow cytometry, for investigating normal and diseased cell states in a wide range of epithelial tissues. Introduction Characterization of protein signaling networks for systems-level analysis of cellular behavior requires the quantification of multiple signaling pathway activities in a multiplex fashion. Previous and current studies of multi-pathway epithelial signaling rely on bulk assays that hinge on the assumption of cell homogeneity in, for example, in vitro cell culture systems. Although useful in revealing coarse-grain biological insights into behaviors exhibited by a majority of cells (Lau et al, 2013, 2012, 2011), these technologies fail to address the complexities exhibited by heterogeneous cell types in vivo. Flow cytometry is a tractable method for detecting and quantifying signal transduction information at single cell resolution (Krutzik et al, 2004; Irish et al, 2004). CyTOF, where the limitation of fluorescence spectral overlap is overcome by the resolution of metal-labeled reagents by mass spectrometry, allows for multiplex sampling of protein signals at a network scale and at single cell resolution (Bendall et al, 2014, 2011). In parallel, newly developed fluorescent dyes and compensation algorithms allow 15-20 parameters to be measured simultaneously using multicolor fluorescent flow cytometry (O’Donnell et al, 2013). A tremendous opportunity for single cell studies lies in expanding quantitative cytometric approaches to epithelial tissues, from which many diseases arise. A significant challenge, however, is the preparation of single cell suspensions from these tissues while maintaining intact cell signaling states. Disruption of epithelial cell junctions during cell detachment perturbs native cell signaling networks (Baum & Georgiou, 2011; Pieters et al, 2012) and can create experimental artifacts that overwhelm native signaling. To date, strategies to quantify epithelial protein signal transduction by cytometry approaches without confounding dissociation artifacts have not been developed. We present a novel method, DISSECT, for preparing single cell suspensions from epithelial tissues for single cell, cytometry-based signaling analyses. We use DISSECT followed by CyTOF to characterize multiple signaling pathway responses in the murine intestinal epithelium following in vivo exposure to TNF-α, a pleiotropic cytokine that plays significant roles in the pathogenesis of Inflammatory Bowel Disease (Colombel et al, 2010), celiac disease (Chaudhary & Ghosh, 2005), and necrotizing enterocolitis (Halpern et al, 2006). In the villus of the duodenum, TNF-α triggers caspase-dependent cell death, creating an epithelial barrier defect that increases exposure of nutrient and microbial antigen to the underlying immune system (Lau et al, 2011; Williams et al, 2013). Remarkably, only a fraction of villus cells undergo apoptosis, and higher levels of cell death cannot be induced by a higher TNF-α dose (Lau et al, 2011). The existence of heterogeneous responses provides a unique opportunity to leverage the natural variation of cells for identifying perturbations that result in desirable cellular outcomes. To decipher heterogeneous responses at single cell resolution, we first provide rigorous, quantitative validation of our single cell approach in comparison to gold-standard lysate-based methods for evaluating both cellular identity and signaling. We then use DISSECTCyTOF to quantify 21 protein and phospho-protein analytes across core signaling pathways at single cell resolution. Quantitative modeling of single cell datasets reveals that a subset of the presumably homogeneous enterocyte population exhibits combinations of signaling responses that confer sensitivity to TNF-α-induced cell death. Our results reveal novel insights into the intricacies of in vivo epithelial cell populations that exhibit significant complexity when perturbed and then observed at single cell resolution. Our approach can be extended to a broad range of complex, heterogeneous epithelial tissues that can be studied via the use of either multiparameter flow cytometry or CyTOF. Results A novel disaggregation procedure for investigating epithelial signaling heterogeneity Tissues in vivo present substantial heterogeneity at the cellular level, as exemplified by the different responses of individual cells to exogenous perturbations. We modeled heterogeneous response in vivo by inducing villus epithelial cell death by systemic TNF-α administration. TNF-α triggered apoptosis only in a third of duodenal villus epithelial cells over a 4-hour time course (Fig. EV1A, B). The remaining cells were not in the process of cell death, as evidenced by the full recovery of intestinal morphology 48 hours after TNF-α exposure (Fig. EV1C). Heterogeneous, TNF-α-induced apoptosis occurred intermittently throughout the length of the villus, and not only at the villus tip as observed in homeostatic cell shedding (Fig. 1A, Fig. EV1D). Furthermore, TNF-α-induced apoptosis appeared to occur solely in a subset of villus enterocytes, as cleaved caspase-3 (CC3) did not co-localize with other epithelial cell type markers (gobletMUC2: Mucin2, tuftDCLK1: doublecortin-like kinase 1, enteroendocrineCHGA: Chromagranin A) (Fig. 1B, Fig. EV1D, E). However, CC3 was co-localized in cells positive for Villin, a protein of enterocyte brush borders, both within the villus epithelium (dying cells) and in the gut lumen (dead cells) (Fig. EV1F). The notion of enterocyte-specific cell death was further supported by increased goblet and tuft cell fractions over time, indicating enrichment of these cell types compared to the remaining enterocytes (Fig. EV1G, H). Although enterocyte cell death occurred heterogeneously in response to TNF-α, the sensing of TNF-α ligand by TNF receptor (TNFR) appeared uniform in these cells. TNFR1 expression was observed on the basolateral membranes of all villus epithelial cells (Fig. 1C, Fig. EV1I), and was reduced in all cells uniformly upon TNF-α stimulation, consistent with internalization of the receptor in direct response to TNF-α binding (Schütze et al, 2008). TNFR2 was expressed at very low levels in the villus epithelium (Fig. EV1I’), supporting previous reports of its minimal role in the villus compartment (Lau et al, 2011). Since TNF-α sensing appeared uniform in all villus epithelial cells, we surmise that heterogeneous TNF-α responses in enterocytes may depend upon differences in signal transduction downstream of receptor binding. A major challenge for exploring signaling heterogeneity in epithelial tissues with cytometry-based methods is the requirement of single cell suspensions. Previous attempts to probe epithelial signaling involved stimulation experiments on single epithelial cells that were already dissociated and outside of their native contexts (Lin et al, 2010). To study single cell signaling in the in situ epithelial context, we first tested whether a single cell disaggregation procedure used routinely for flow sorting epithelial cells (Magness et al, 2013) (which we referred to as “the conventional method”) preserves native signaling in single cell suspensions. Briefly, the intestinal epithelium was mechanically retrieved after the intestine was acquired, washed, and longitudinally opened. The epithelium was then digested enzymatically (~10 min) and then filtered into a single cell suspension. A standard fix-perm procedure for phospho-flow was then performed, followed by phospho-specific antibody staining and cytometry analysis (Krutzik et al, 2004). Quantitative immunoblotting analysis on fresh intestinal tissue lysates was used as positive control. A head-to-head assessment using the same antibodies was performed by comparing median intensities from single cell flow cytometric data to integrated intensities of bands from immunoblots, which reflect cell averages in tissue lysates. This comparison demonstrated that signal transduction induced by TNF-α was not maintained with the conventional disaggregation method, as assessed by both early (p-ERK1/2, p-C-JUN) and late (p-STAT3) signals (Fig. EV2). A previous study suggested that signaling perturbations from tryptic disaggregation can be eliminated by performing digestion in live cells at low temperatures (Abrahamsen & Lorens, 2013). We tested the effect of enzymatic digestion by performing lowyield single cell disaggregation on live tissues (Appendix Fig. S1A), using gentle mechanical dissociation without any enzymes; however, signal transduction was still not preserved (Appendix Fig. S1B). Disaggregation of an intact epithelium into single cells perturbs components of epithelial cell junctions that play many roles in signaling modulation. Such disruption in live tissue may dynamically alter signaling pathways and produce experimental artifacts. To adapt single cell signaling analysis for epithelial tissues, we developed DISSECT, a single cell dissociation method that preserves intact signaling. After the epithelium was retrieved from the animal, it was immediately fixed to maintain cellular signaling states. The epithelium was then subjected to acetone permeabilization and antigen retrieval by a detergent solution, followed by staining and an additional fixation step to crosslink antibodies onto their epitopes. Stained epithelium was then disaggregated into single cells enzymatically followed by gentle mechanical agitation (Fig. 1D). Retrieval of single cells and their yields were robustly verified, with cells prepared by DISSECT retaining a native columnar morphology, versus the round morphology arising from the conventional method (Appendix Fig. S2, Fig. EV3A, B). Specifically, quantitative yields of single cells from DISSECT were higher than those from the conventional approach, where cell clumping induced by methanol and pronounced adhesion of single cells to plastic ware resulted in cell loss (Fig. EV3C). We tested whether native signaling is maintained throughout the DISSECT process, again by direct comparison with “gold standard” approaches performed on the same tissues. Activation of p-C-JUN and p-STAT3 were detected at 0.5 hour and 4 hours, respectively, in single cells by immunofluorescence microscopy, mirroring intact tissue staining (Fig. 2A). Singleplex flow cytometry on prepared single cell suspensions enabled the quantification of signal transduction at single cell resolution, whose median values from single cell distributions can be compared to lysate-based quantitation (Fig. 2A’). We observed upregulation of p-C-JUN early (0.5 hour) and p-STAT3 late (2 hour) at the population level, matching previously observed dynamics of these two TNF-αactivated pathways (Lau et al, 2011). Furthermore, median data derived from single cells prepared using DISSECT over multiple replicates qualitatively matched immunoblotting data from lysates prepared from the same tissue (Fig. EV2), in stark contrast to single cells prepared using the conventional method. Furthermore, preservation of signals using DISSECT was not furthered improved by perfusing the animal beforehand with fixative, indicating that our method of tissue collection does not significantly perturb native signaling (Appendix Fig. S3). By verifying the performance of DISSECT in technical and biological replicates (Appendix Fig. S4), we conclude that our procedure is robust for maintaining native signaling during single cell disaggregation. DISSECT allows phenotypic cell profiling of complex epithelial tissues A potential limitation of DISSECT is the possible degradation of proteins at the cell surface, thus limiting our ability to identify cell types using cell surface markers. To ensure that the DISSECT approach can preserve cell surface antigen staining for cell-type identification, we evaluated canonical markers for leukocytes and other epithelial cell types using flow cytometry in our single cell preparations. CD45+ cells in the intestinal lamina propria can be readily detected and increased 2 hours after TNF-α stimulation, similar to what we observed previously for immune cell types (Appendix Fig. S5A) (Lau et al, 2012). Specifically, we detected different populations of villus epithelial cells using goblet (CLCA1 Calcium-activated chloride channel regulator 1), enteroendocrine (CHGA), and tuft (DCLK1) cell markers (Fig. 2B). The proportion of differentiated cells detected in the intestine matched previous reports, with goblet cells at ~10% and increasing from the duodenum to the ileum (Rojanapo et al, 1980; Paulus et al, 1993; Imajo et al, 2014; van der Flier & Clevers, 2009; Wright & Alison, 1984), enteroendocrine cells at ~1% (Gunawardene et al, 2011; Cheng & Leblond, 1974), and tuft cells at ~1% (Gerbe et al, 2012). Imaging-based quantification of the same tissues also confirmed these results (Fig. 2B’, B’’). We further tested whether our method can detect crypt stem cells using the cell surface marker LRIG1 (Appendix Fig. S5B) (Powell et al, 2012). Isolation of colonic crypts followed by DISSECT and flow cytometry allowed for the identification and quantification of crypt base cells, which segregate away from Na/ATPase+ differentiated cells (Fatehullah et al, 2013) (Appendix Fig. S5C). The proportion of LRIG1+ cells matched what was previously reported (~30% in the colonic crypt) using the same antibody (Poulin et al, 2014). In addition, TNF-α-induced signaling can be detected in single cells isolated from colonic crypts (Appendix Fig. S5D), as well as from colonic tumors (Appendix Fig. S5E) using flow cytometry following DISSECT. To test the general applicability of DISSECT in other epithelial tissues, we induced proliferation in the collecting ducts of the kidney and hepatocytes in the liver using an unilateral ureteral obstruction (UUO) model and a partial hepatectomy model, respectively. GFP+ cells from the Hoxb7cre;mT/mG mouse labels collecting duct cells, which can be identified by flow cytometry postDISSECT (Fig. EV4A, A’). UUO-induced injury triggered proliferative responses to various degrees in different mice, which correlated with p-RB proliferative signaling (Fig. EV4B) (Giacinti & Giordano, 2006). Furthermore, after partial hepatectomy, BrdU-labelled hepatocytes (Fig. EV4C, D) were enriched for p-RB signaling during the recovery phase (Fig. EV4E). These results demonstrate DISSECT to be a valid, reliable approach for disaggregating a variety of heterogeneous epithelial tissues into single cell suspensions for cytometry-based signaling analysis. DISSECT preserves signal transduction across a wide range of signaling pathways in epithelial tissues We expect comparable quantitative approaches to have relatively comparable signal-tonoise detection. With regards to noise, we compared the standard deviation of signals generated from biological replicates using different quantitative approaches. Results generated by DISSECT followed by flow cytometry matched with those obtained by lysate-based ELISA assays and quantitative immunofluorescence imaging, demonstrating that these assays pick up comparable levels of noise (Appendix Fig. S6). With regards to signal, we performed rigorous, quantitative comparisons of TNF-α-induced signaling measurements between DISSECT-flow cytometry and two gold standard methods: quantitative immunofluorescence imaging (Fig. 3) and quantitative immunoblotting (Appendix Fig. S7). A summary of how we derived quantitative information from each of the three methods is documented in Appendix Figure S8. The same set of antibodies was used for all three methods to evaluate protein states, such as phosphorylation and cleavage, that act as direct surrogates of signaling pathway activation. Three cohorts of mice (30 samples) were used for each analysis, and tissues from each animal were split three ways for different types of analyses. Because lysate-based approaches assess the average of all cell types in whole tissue, our cytometry analyses were also performed in a bulk cell population manner to enable direct comparison between approaches. To sample a wide dynamic range, we leveraged tissues from the duodenum and ileum (which exhibit differential TNF-α signaling responses), as well as from different time points post TNF-α exposure to generate quantitative correlation analyses. Ten out of eleven protein analytes generated statistically significant correlations between DISSECT-flow quantification and imaging quantification (6 out of 6 with quantitative immunoblotting) (Fig. 3, Appendix Fig. S7). Combined correlation analyses using all protein analytes resulted in a highly significant correlation (p<0.0001) between DISSECT-flow and imaging data, and between DISSECT-flow and immunoblotting data. The Pearson’s coefficients of comparing DISSECT-flow to imaging and immunoblotting were 0.72 and 0.81, respectively. Factors that contribute to the imperfect correlation include inherent experimental noise and differences in quantification between each of the methods, which will be discussed below. Furthermore, for a truly unbiased analysis, we did not exclude obvious data outliers that affected the normalization procedure, which can skew relatively small datasets and can subsequently weaken the correlations. Nevertheless, our conservative approach for validation still generated highly significant (p<0.0001) correlations. These results demonstrate the validity of DISSECT to preserve native signaling during single cell disaggregation, and to generate single cell level data, when aggregated as populations, detect similar signal-to-noise as “gold standard” population-based methods. DISSECT application of CyTOF identifies a differentially signaling enterocyte subpopulation that is sensitized to TNF-α-induced cell death A 21-analyte CyTOF panel of heavy-metal labeled reagents specific for epithelial signaling was generated (Appendix Table S1). Twenty one-plex CyTOF analysis was performed on three cohorts of mice subjected to a time course of acute TNF-α exposure, giving rise to average early and late signaling results that matched with flow cytometry, imaging and quantitative immunoblotting (Fig. 4A). We used single cell CyTOF data to first reaffirm TNF-αinduction of cell death strictly within the duodenal enterocyte population. Indeed, CC3 did not colocalize with other epithelial cell type-specific markers (CK18: cytokeratin 18 – secretory subset, CLCA1 goblet, CHGAenteroendocrine, CD45 leukocytes) (Fig. 4B, C compared to Fig. EV1E). The few double positive cells are not cell clusters (Appendix Fig. S9). The fraction of differentiated cell types detected again matched published results (Rojanapo et al, 1980; Paulus et al, 1993; Imajo et al, 2014; van der Flier & Clevers, 2009; Wright & Alison, 1984; Gunawardene et al, 2011; Cheng & Leblond, 1974; Gerbe et al, 2011), as well as flow and imaging data we obtained previously (Fig. 2B, 4B). To identify subpopulations of enterocytes with distinct signaling activities indicative of cell death, we used t-SNE (t-Distributed Stochastic Neighbor Embedding) to visualize multiplex single cell data in 2 dimensions while maintaining dissimilarities between cells in multidimensional data space (Fig. 4D, Dataset EV1) (Amir et al, 2013). We again focused on the one-hour time point to characterize active signaling cells undergoing cell death. t-SNE analysis allowed groupings of different functional cell types based on combinations of signaling and cell identity markers. In addition, a distinct population of CC3+ enterocytes was identified. We used manual gating on t-SNE space to supervise a partial least squares discriminant (PLSDA) model to categorize enterocytes undergoing cell death against living enterocytes. Classification based upon calibration signaling data in 2-latent variable PLSDA space to predict CC3 expression resulted in an area (AUC) of 0.92 under the receiver of operating characteristic (ROC) curve, indicative of high sensitivity and specificity (Fig. 4E). We then cross-validated our model by repeatedly withholding 10% of the data using random, venetian blind, and block selection. Our cross-validation model yielded similar prediction power (ROC AUC = 0.92) compared to our calibration model due to the high number of data points used for fitting a model with a relatively limited set of parameters, which dramatically lowers the prospects of overfitting. We used the discriminant coefficients (β) of our PLSDA model to select signaling features that were informative for classification. Using 10000-fold permutation testing, we generated β distributions around zero and determined the probability for obtaining our model coefficients. The four coefficients with the lowest p–values were p-P38, p-CREB, p-ERK, and CK20 (Fig. 4F). Another method for feature selection using Variable Importance in Projection (VIP) scores also identified the same four variables (Fig. 4G). We overlaid these four variables onto t-SNE plots to determine their ability to predict CC3 expression (Fig. 4H). While individual variables positively or negatively correlated with the CC3+ population, they were incapable of clearly discerning this population from other cellular populations (Fig. 4I). Linearly combining these four variables without scaling allowed for clear identification of CC3+ enterocytes (Fig. 4J), indicating that combinatory activities of multiple signaling pathways contribute to a “signaling code” that implicate cell death. More importantly, the same experimental and computational analysis applied to three different cohorts of mice selected the same set of four variables that identify CC3+ enterocytes (Fig. 5, Dataset EV2, EV3). In addition, other β coefficients besides the top four variables also followed the same trend of positive or negative correlation with CC3 in different mouse cohorts. These results indicate that DISSECT followed by CyTOF is a highly reproducible method to accurately characterize single cell behavior using multi-pathway signaling parameters. Divergently responding enterocytes are neighbors within the intestinal epithelium Having a signaling fingerprint that classifies dying and non-dying enterocytes allows us to identify divergent signaling mechanisms that significantly affect intestinal physiology. Specifically, we chose to investigate divergent p-ERK signaling in the intestinal epithelium, which occurred in the surviving, but not in the dying, cell population. p-ERK activation in surviving enterocytes was also heterogeneous, which prompted us to envision spatial patterns of p-ERK activity that conferred survival. Whole-mount imaging of whole villus at 1 hour post TNF-α exposure revealed a “flower petal” ring-like pattern of epithelial p-ERK signaling, with five or six p-ERK positive cells surrounding a p-ERK negative area (Fig. 6 A, Fig. EV5A, yellow arrows). Co-staining with CC3 revealed that in many cases, the dying CC3+ cells occupied the central area surrounded by p-ERK+ neighbors (Fig. 6B, Fig. EV5B, yellow arrows). In other cases, the dying CC3+ cell has already been extruded from the epithelium, leaving an apoptotic rosette surrounded by p-ERK+ cells undergoing contraction-dependent closure (red arrow). Furthermore, the ratios of CC3+ dying cells and p-ERK+ enterocytes in 3 cohorts of mice were 1:4.56, 1:6.04, 1:4.73, respectively, supporting that the immediate neighbors of the dying cell activated p-ERK signaling (Fig. EV5C, D). Imaging of tissue sections also corroborated that dying cells were flanked by p-ERK+ cells (Fig. EV5E), although the phenomenon was harder to visualize in 2 dimensions. We surmise that the dying cell signals to neighboring cells nonautonomously to activate a cell survival program, in order to prevent large swaths of contiguous epithelium from dying and to prevent unrecoverable barrier defects. Thus, we tested the effect inhibiting p-ERK signaling using the allosteric MEK inhibitor PD0325901 (Fig. EV5F). Inhibition of p-ERK signaling affected the latency of the cell survival program such that epithelial apoptosis occurred immediately following TNF-α exposure, which resulted in a higher number of dying cells in total (Fig. 6C). Inhibition of P38 alone minimally affected TNF-α induced apoptosis (Fig. EV5G), but was able to partially normalize early apoptosis due to MEK inhibition (Fig. 6C), consistent with P38’s pro-apoptotic role that is context-dependent. To our knowledge, this is the first reported observation of this “flower petal” pattern of p-ERK activation in response to TNF-αinduced cell death in epithelial tissue. This new finding demonstrates the applicability of our single cell signaling experimental platform, in conjunction with data analysis, to reveal novel, non-cell autonomous responses in complex heterogeneous epithelia.